The dataset contains the following variables:
name: the name of the meteorite (typically a location, often modified with a number, year, composition, etc) id: a unique identifier for the meteorite
nametype: one of:
recclass: the class of the meteorite; one of a large number of classes based on physical, chemical, and other characteristics (see the Wikipedia article on meteorite classification for a primer)
mass: the mass of the meteorite, in grams
fall: whether the meteorite was seen falling, or was discovered after its impact; one of:
year: the year the meteorite fell, or the year it was found (depending on the value of fell)
reclat: the latitude of the meteorite’s landing
reclong: the longitude of the meteorite’s landing
GeoLocation: a parentheses-enclose, comma-separated tuple that combines reclat and reclong
We’ll clean out incorrect reccomended lat/long coordinates, we’ll also filter for USA via using a bounding box for lat/long within the most north/south/east/west bounds of the continental US.
meteorites <- meteorites %>%
dplyr::filter(year >= 860 & year <= 2016 ) %>% # filter out weird years
dplyr::filter(reclat != 0 | reclong != 0) %>%
dplyr::arrange(desc(year)) %>%
dplyr::filter(!is.na(mass)) %>%
dplyr::filter(mass > 0) %>%
dplyr::mutate(year = as.numeric(year))
# src for bbox
# https://gist.github.com/graydon/11198540
meteorites_clean <- meteorites %>%
dplyr::filter(reclat <= 49.3457868 & reclat >= 24.7433195) %>% #filter for north and south
dplyr::filter(reclong >= -124.7844079 & reclong <= -66.9513812) # filter for east and west
pal <- colorNumeric(palette = "RdYlBu", domain = meteorites_clean$mass, reverse = T)
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal(mass),
radius = 3,
stroke = FALSE, fillOpacity = 0.5) %>%
addLegend("bottomright", pal = pal, values = ~mass,
title = "Mass (grams)",
labFormat = labelFormat(suffix = " g"),
opacity = 1
)
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addHeatmap(lng = ~reclong, lat = ~reclat,
gradient = "viridis",
minOpacity = 0.1,
cellSize = 1,
radius = 3, blur = 15)
Here we’ll visualize the meteorites via whether they were spotted on the ground after impact (“found”) or while in the mid-fall ("fell)
Overwhelmingly we can see most meteorites were spotted while on the ground after-impact was made, we have modified the size argument to help make the quantity more acute on our map.
pal <- colorFactor(c("#6EC5E9", "#FF1C51"), domain = c("Fell", "Found"))
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal(fall),
radius = ~ifelse(fall == "Found", 1, 3),
stroke = FALSE, fillOpacity = 1) %>%
addLegend("bottomright", pal = pal, values = ~fall,
title = "Found Status",
opacity = 1)
Distribution of Meteorites by Years
# Create a continuous palette function
pal <- colorBin(
palette = "Spectral",bins = 5,
domain = meteorites_clean$year, reverse = T)
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal (year),
radius = 3,
stroke = FALSE, fillOpacity = 1) %>%
addLegend("bottomleft",
colors =c("#2B83BA", "#ABDDA4", "#FFFFBF","#FDAE61","#D7191C"),
labels = c("1600 - 1700", "1700 - 1800", "1800 - 1900", "1900 - 2000", "2000 - 2016"),
title = "Year Found",
opacity = 1)
Most of the meteorites were discovered within last century (1900 - 2000)
Frequency distribtion by year
hchart(meteorites_clean$year,
color = "Purple",
name = "year") %>%
hc_yAxis(title = list(text = "Frequency (n)")) %>%
hc_xAxis(title = list(text = "<b> Year Observed </b>")) %>%
hc_add_theme(hc_theme_monokai())
We’ll divide the observation into a grid (h = 5, width = 10) and count the points that lie within each quadrat
# get the window for the data, i.e. the boundaries
W <- owin(xrange = c(min(meteorites_clean$reclong),
max(meteorites_clean$reclong)),
yrange = c(min(meteorites_clean$reclat),
max(meteorites_clean$reclat)))
# select just the lat/long columns
pp1 <- meteorites_clean %>%
dplyr::select(reclong,
reclat)
# convert to ppp object
meteorites_ppp <- as.ppp(pp1, W = W)
Q <- quadratcount(meteorites_ppp, nx= 5, ny=10)
# plot results
plot(meteorites_ppp, pch=20, cols="#FF1C51", main=NULL) # Plot points, with the tallies
plot(Q, add=TRUE) # Add quadrat grid
meteorites_ppp.km <- rescale(meteorites_ppp, 1000, "km")
# Compute the density for each quadrat (in counts per km2)
Q <- quadratcount(meteorites_ppp.km, nx= 6, ny=3)
Q.d <- intensity(Q)
# Plot the density
plot(intensity(Q, image=TRUE), main=NULL, las=1) # Plot density raster
plot(meteorites_ppp.km, pch=20, cex=0.6, col="#39FF14", add=TRUE) # Add points
K1 <- density(meteorites_ppp.km) # Using the default bandwidth
plot(K1, main=NULL, las=1)
contour(K1, add=TRUE)
We’ll apply three difference measures for distance bases analysis:
Computing the ANN distances between imapct sites. Compute the average first nearest neighbor distance (km), set k = 1
# set k - 1 for the first nearest neighbor distance
mean(nndist(meteorites_ppp.km, k = 1))
## [1] 0.0002431465
Compute the average second nearest neighbor, set k = 2
mean(nndist(meteorites_ppp.km, k = 2))
## [1] 0.0003829143
ANN <- apply(nndist(meteorites_ppp.km, k=1:2000),2,FUN=mean)
# generate an ANN vs neighbor order plot
plot(ANN ~ eval(1:2000), type="b", main="ANN Plot", las=1)
The bottom axis shows the neighbor order number and the left axis shows the average distance (km). We can see a linear trend via
Computing the K function
K <- Kest(meteorites_ppp.km)
plot(K, main=NULL, las=1)
This plot is visualizing the different estiamtes of \(K\) depending on the edge correction. The edge corrections plotted are listed in the legend.
\(K_{iso}\) refers to an isotropic implementation
\(K_{trans}\) refers to an translate implementation
\(K_{bord}\) refers to an border implementation
\(K_{pois}\) refers to an border implementation
Here we’ll estimate the meteorite point process intensity using the mass
In the following example, a Starbucks store point process’ intensity is estimated following the population density raster covariate
Running an average nearest neighbor analysis for impact sites, in this we assume a uniform point density across space
ann.p <- mean(nndist(meteorites_ppp.km, k =1))
ann.p
## [1] 0.0002431465
The observed average nearest neighbor distance is 0.0002431465 km (0.24 m).
Generate the distribution of an expected
Useful good here
xy <- meteorites_clean[,c(9,8)]
meteorites_clean <- SpatialPointsDataFrame(coords = xy, data = meteorites_clean,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
print(nni(meteorites_clean))
## Warning: data contain duplicated points
## $NNI
## [1] 0.6192556
##
## $z.score
## [1] -30.13813
##
## $p
## [1] 1.534619e-199
##
## $expected.mean.distance
## [1] 0.3926433
##
## $observed.mean.distance
## [1] 0.2431465
Our NNI score is less than 1, we can consider this as a good indication of clustering and can reject the Null Hypothesis that the points are uniform in pattern
In this section we’ll tally the number of meteorites via the State in which they land to gather some areal statistics.
# Downloading US State polygons
states <- states(cb = TRUE)
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
# Convert to sf data
points <- st_as_sf(meteorites_clean)
# Set CRS
states <- st_transform(x = states, crs = 4326)
points <- st_transform(x = points, crs = 4326)
# Intersection between polygon and points ---------------------------------
intersection <- st_intersection(x = states, y = points)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Count number ofintersections per polygon (i.e. states)
states$pt_count <- lengths(st_intersects(states, points))
# Interactive Map
pal <- colorNumeric(
palette = "viridis",domain = states$pt_count, reverse = F, alpha = T)
labels_info = sprintf(
"<h2> %s </h2>
<hr>
Meteorites: %s ",
states$NAME,
states$pt_count
) %>% lapply(htmltools::HTML)
leaflet(states) %>%
addProviderTiles("CartoDB.DarkMatter") %>%
addPolygons(fillColor = ~pal(pt_count),
weight = 2,
opacity = 0.25,
color = "white",
dashArray = "1",
fillOpacity = 0.6,
highlight = highlightOptions(
weight = 3,
color = "#8ffcff",
fillColor = "#8ffcff",
dashArray = "1",
fillOpacity = 0.6,
bringToFront = TRUE),
label = labels_info,
labelOptions = labelOptions(
style = list("font-family" = "arial",
"font-style" = "bold",
"color" = "#FF1C51",
"background" = "black",
"border-color" = "#FFD300")
)
) %>%
setView(-98.5795, 39.8282, zoom=3) %>%
addLegend("bottomright",
pal = pal,
values = ~pt_count,
title = "Meteorites <hr>",
opacity = 1)